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ABSTRACT 

We describe a new method to extract spectra of stars from observations of crowded stellar fields with integral field spectroscopy 
(IFS). Our approach extends the well-established concept of crowded field photometry in images into the domain of 3-dimensional 
spectroscopic datacubes. The main features of our algorithm are: (1) We assume that a high-fidelity input source catalogue already 
exists, e.g. from HST data, and that it is not needed to perform sophisticated source detection in the IFS data. (2) Source positions 
and properties of the point spread function (PSF) vary smoothly between spectral layers of the datacube, and these variations can be 
described by simple fitting functions. (3) The shape of the PSF can be adequately described by an analytical function. Even without 
isolated PSF calibrator stars we can therefore estimate the PSF by a model fit to the full ensemble of stars visible within the field of 
view. (4) By using sparse matrices to describe the sources, the problem of extracting the spectra of many stars simultaneously becomes 
computationally tractable. We present extensive performance and validation tests of our algorithm using realistic simulated datacubes 
that closely reproduce actual IFS observations of the central regions of Galactic globular clusters. We investigate the quality of the 
extracted spectra under the effects of crowding with respect to the resulting signal-to-noise ratios (S/N) and any possible changes in 
the continuum level, as well as with respect to absorption line spectral parameters, radial velocities and equivalent widths. The main 
effect of blending between two nearby stars is a decrease in the S/N in their spectra. The effect increases with the crowding in the field 
in a way that the maximum number of stars with useful spectra is always ~0.2 per spatial resolution element. This balance breaks 
down when exceeding a total source density of ~1 significantly detected star per resolution element. We also explore the effects of 
PSF mismatch and other systematics. We close with an outlook by applying our method to a simulated globular cluster observation 
with the upcoming MUSE instrument at the ESO-VLT. 

Key words. Methods: data analysis - Techniques: imaging spectroscopy - Galaxy: globular clusters 
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1. Introduction 

The observation of spatially resolved, but densely populated stel- 
lar fields such as star clusters or the inner regions of very nearby 
galaxies poses a severe challenge because of crowding. At finite 
angular resolution each star is represented by the point spread 
function (PSF) characteristic of the system. Observable is only 
the superposition of many scaled and shifted PSFs, according to 
brightness and position of the stars in the field of view. It is then 
a major task to disentangle the various overlapping PSFs in or- 
der to measure the corresponding stellar fluxes as accurately as 
possible. 

In imaging photometry, this challenge was addressed already 
many years ago, and dedicated analysis tools have been devel- 
oped to perform what is commonly termed 'crowded field pho- 
tometry', daophot (Stetson 1987) is the most widely used of 
such tools and has been applied extensively and with great suc- 



* Based on observations collected at the Centro Astronomico His- 
pano Aleman (CAHA) at Calar Alto, operated jointly by the Max- 
Planck Institut fur Astronomie and the Instituto de Astroffsica de An- 
dalucia (CSIC) 

** Based on observations made with the NASA/ESA Hubble Space 
Telescope, and obtained from the Hubble Legacy Archive, which 
is a collaboration between the Space Telescope Science Institute 
(STScI/NASA), the Space Telescope European Coordinating Fa- 
cility (ST-ECF/ESA) and the Canadian Astronomy Data Centre 
(CADC/NRC/CSA). 



cess. At the core of daophot (and similar software) is a code 
that fits a PSF model simultaneously or in succession to all stars 
within a given field. The resulting photometry is unbiased with 
respect to overlapping stellar images (within the accuracy of the 
PSF model and as long as the list of fitted stars is complete), so 
that even heavily blended stars can be recovered. 

However, the astrophysically desirable move from photome- 
try to spectroscopy is still fraught with difficulties when it comes 
to crowded stellar fields. Traditional long-slit (or multi-slit) 
spectroscopy is clearly limited by the amount of blending within 
the slit aperture, which under conditions of heavy crowding may 
be almost impossible to control. The same is true for fibre-fed 
multiplexing spectroscopy using single apertures. Consequently, 
investigations using this type of equipment have typically re- 
stricted themselves to regions with modest crowding, such as 
the outer parts of star clusters or galaxies, and/or focused on the 
brightest sourc es in the fields o f interest. A nota ble excep tion 
is the study by |van der Marel et aT] ( |2002| ) and |Gerssen et al. 
( 2003) who performed HST long-slit spectroscopy of the center 
of the globular cluster Ml 5. Here the improved spatial resolu- 
tion of HST helped in dealing with the limitations of traditional 
spectroscopy. 

Potentially much more powerful in this domain is the direct 
combination of imaging and spectroscopy. 'Integral Field Spec- 
troscopy' (IFS) - often also called 3D or IFU spectroscopy (IFU 
= integral field unit) - has matured into a widely used observ- 
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ing technique over the last decade. It is particularly powerful for 
extended objects such as galaxies, providing access to spectro- 
scopic information across the full spatial extent of the target. In- 
tegral field spectroscopy has strongly advance d our knowledge 
of galaxies, both in the local universe (e.g., Cappellari et al. 



201 1 ; de Zeeuw et al. 2002 ) and at higher redshifts (e.g., |Forster 
Schreiber et al. 12009) ). 



Beyond the ability to trace the variation of observable prop- 
erties across the field of view (FoV), IFS data also provide the 
user with an unprecedented access to single out individual ob- 
jects surrounded by (and possibly blended with) other sources. 
This can be done in increasing levels of sophistication. A 
straightforward approach in stellar fields is to identify spatial 
pixels that are dominated by the contribution of a single star 
and obtain a spectrum of that star by summing the spectra in 
only those pixels. An example is the recent work by |Evans et al.| 
( |2011| l, who used IFS to obtain spectra of massive stars in the 
Tarantula nebula (30 Doradus). A similar treatment was applied 
to n ear-infrared IFS observ atio ns of stars in the G alactic centre 
by |Eisenhauer et al.| ( |2005) and |Pfuhl et aTlpOTT] ). While scien- 
tifically very successful, this approach is obviously still limited 
to moderate crowding conditions. 

In order to facilitate the extraction of the spectra, knowl- 
edge of the spatial shape of the sources must be incorporated. 
Especially when dealing with barely resolved or yet unresolved 
sources, a proper knowledge of the PSF becomes extremely im- 
portant. Unfortunately, the fields of view of existing IFS in- 
struments are so small that they often do not contain any undis- 
turbed PSF calibrator. In order to proceed it thus becomes neces- 
sary to estimate the PSF directly from the data. Several authors 
have presented solutions of this task, although usually restricted 
to specific science applications. For the severely background- 
limited spectrophotometry of planetary nebulae (PN) in the 
bulge of M3 l, |Rothetal.lp004| ) estimated the PSF from a bright 
PN image in [O in] and fitted this to each wavelength bin of the 
datacu be. |Fabrika et al.| ( |2005| ) applied the cplucy algorithm 
(Lucy 1994) with input from high spatial resolution HST im- 
ages to extract spectra of LBV candidates in M33, while fully 
accounting for a highly variable unresolved background. Tech- 
nical details of these two applications are discussed in Becker 
et al.| ( |200"4) . In a similar spirit, the bright nuclear emission lines 
of quasars can be used to 'self-calibrate' the PSF for the analy- 
sis of quasar host galaxies (Christen sen et al.||2006[ |Husemann| 
^TaTl|MTll |Jahnke et al.||2004) . However, due to the lack of 
emission line point sources, these procedures cannot be applied 
to arbitrary stellar fields. Another limitation lies in the fact that 
emission lines provide the PSF only at specific wavelengths, and 
an extrapolation to the full spectral range may be doubtful. 

A more general approach was explored by Wisotzki et al.| 
2003) who deblended the four overlapping components of the 



gravitationally lensed QSO HE 0435-1223 from IFS data by as- 
suming a purely analytic (but wavelength-dependent) PSF shape 
that was iteratively optimised. The main limitation of this par- 
ticular solution was its computational inefficiency, which would 
make an application to fields with many more sources pro- 
hibitive. Sanch ez et al.| ( |2007) extended this approach to non- 
point sources in an IFS study of the galaxy cluster Abell 2218, 
where they used the morphological information on individual 
galaxies obtained from high-resolution imaging to deblend the 
overlapping data into individual galaxy spectra. 

In this paper we present a new algorithm for the spectropho- 
tometric analysis of generic crowded stellar fields observed by 
integral field spectroscopy. The algorithm determines a fully 
self-calibrated wavelength-dependent PSF model which is sub- 



sequently fitted to the entire datacube. In many aspects the al- 
gorithm is an extension of the well-established 'crowded field 
photometry' approach, but the spectroscopic dimension requires 
the addition of some genuinely new features. Note that while we 
were completing this paper, [Soto et al.| ( [2012) published an anal- 
ysis of IFS observations in the Galactic bulge which includes 
PSF estimation and bears some generic resemblance to our ap- 
proach. Unfortunately, not many details are provided in that pa- 
per on the algorithm itself and its performance. 

The purposes of the present paper are twofold. Firstly, the 
paper provides the methodical foundations for a number of sub- 
sequent articles focusing on the central regions of globular clus- 
ters (Kamann et al., in prep.). Here we use some of those data for 
illustration purposes only. Secondly, we also want to highlight 
the huge potential of IFS for crowded field spectroscopy in gen- 
eral, especially in view of the upcoming wide-field panoramic 
integral-field spectrograph MUSE at the ESO-VLT. 

The paper is organized as follows. In Sect. [2] we start with an 
exposition of the basic considerations that went into the develop- 
ment of our method. An overview of the observational data that 
motivated this study, and of the simulated data that we created 
to test our algorithm, is given in Sect. [3] The analysis algorithm 
itself is presented in detail in Sect. [4] Section [5] describes the ex- 
tensive tests that we have carried out to validate its performance. 
A discussion of potential sources that could systematically influ- 
ence the performance of the algorithm is presented in Sect. [6] In 
Sect. [7] we demonstrate briefly the application of our method to 
data that will soon be obtained with the MUSE instrument. We 
wrap up our conclusions in Sect. [8] 

2. Basic concepts 

In a simplified picture, an integral field datacube can be consid- 
ered a sequence of monochromatic images, hereafter called lay- 
ers. A straightforward approach would then be to use existing 
methods to analyze the data cube layer by layer, i.e. to perform 
crowded field photometry individually on each layer, and com- 
bine the results afterwards. Yet, as we shall demonstrate in the 
following, such an approach does not use the full potential of 
IFS data. 

The development of our new algorithm was guided by the 
following thoughts: 

(1) For many objects, like most globular clusters or several 
nearby galaxies, high resolution images are already available 
(e.g., the HST/ ACS survey of galactic Globular Clusters by |Sara- 
jedini et al. 2007). The depth reached by these images is usu- 
ally sufficient to assume that all stars that can be resolved with 
a (seeing-limited) IFU have been detected, and that their relative 
positions are known. Thus, an inventory of stars in the observed 
field already exists, and there is no need to perform sophisticated 
source detection on the integral field data. But of course the 
question remains which sources from an available catalogue can 
be recovered in the datacube. Given the typically lower spatial 
resolution of IFS data, only a subset of the catalogued sources 
will be accessible to the analysis. We describe in Sect. |4.5| how 
we derive an optimal subset of sources from a catalogue that can 
be extracted reliably. 

(2) The individual monochromatic layers of a data cube are 
not independent from one another. Because all of them were ob- 
served simultaneously, temporal effects like seeing variability, 
atmospheric dispersion and instrument flexure affect all images 
in a way that properties such as the PSF or source positions are 
interconnected. While they may vary from layer to layer, these 
variations will be smooth with wavelength. We thus should use 
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the whole datacube to determine one (wavelength dependent) 
PSF model and one set of (wavelength dependent) coordinates. 
The potential of IFS to obtain coordinates with a precision be- 
yond that achievable with a single image at the same spatial res- 



olution has been demonstrated previously (e.g., Mediavilla et al. 
[T998l|Wisotzki et al.|2003| ). 

(3) For the current generation of integral field spectrographs, 
the relatively small field of view often implies that sufficiently 
bright and isolated PSF calibrator stars are not available within 
a science exposure. Dedicated (i.e. non- simultaneous) observa- 
tions of isolated stars are of very limited use for accurate PSF 
modelling because of the temporal variations in the PSF. In the 
future it may become possible to pre dict the PSF from ad aptive 
optics wavefront sensor data (e.g., Jolissaint et al. 2010]). For 
the present, however, we are bound to estimate the PSF directly 
from the crowded field science data. Nevertheless, this informa- 
tion is actually there - not in individual stars, but encased in the 
ensemble of all stars in the field. In Sect.|4.2|we present an itera- 



tive procedure to construct a global PSF model by simultaneous 
fitting of the full ensemble. 

(4) Besides the spectra of individually resolved stars, an IFS 
data cube may also contain a background of fainter and possi- 
bly unresolved stars, and/or nebular emission. This quasi-diffuse 
component can be a major nuisance and source of systematic er- 
rors during the extraction of stellar spectra, but it is also possible 
that the scientific focus is mainly directed towards this compo- 
nent. For these reasons we made a special effort to find an ap- 
propriate solution to account for such a non-trivial background 
component in the analysis. This is explained in detail in Sect [4] 

3. Data 

3.1. Observations of globular clusters 

We recently collected 3D spectroscopy data of a sampl e of 
Galactic Globular Clusters with the instruments PMAS (Roth 
eTaLl|2QQ5] ) and ARGUS ( [Pasquini et aLl[20Q2| ). In each clus- 
ter we mapped the central region out to a distance comparable 
to the core radius of the cluster. The main scientific aims are 
(1) a search for binary stars, and (2) to place constraints on the 
possible presence of intermediate-mass black holes in the clus- 
ters. The results from these observations will be presented in 
forthcoming papers (Kamann et al, in prep.). 

The instrument setups in these observations varied a little 
from campaign to campaign, but were generally chosen to facili- 
tate a kinematic analysis given the expected velocity dispersions 
of the clusters. The typical spectral resolution R - XI AX was 
-8000, and the spectral range was targeted at the Calcium triplet 
0U8498 A, 8542 A, 8662 A). In order to sample the PSF in 
the best possible way we always selected the smallest available 
spaxel scale (spaxel = spatial pixel), 0Y3 for ARGUS and 075 for 
PMAS. With 16 x 16 spaxels for PMAS and 22 x 14 for ARGUS 
we thus cover an area per pointing of 8" x 8" and 676 x 472, 
respectively. The seeing conditions were variable, with ~077 in 
the best and 174 in the worst case. Some whitelight images are 
presented Fig. [T] to give an impression of the data. The crucial 
dependence of the data quality on the seeing during the observa- 
tions is clearly visible. 

3.2. Simulated data 

To assess the the performance of our algorithm and to validate 
the quality of the deblended spectra we carried out extensive 
simulations. In line with our main scientific interest we focus 



on Globular Clusters, but the results are of course applicable to 
any sort of crowded stellar fields. At any rate the central regions 
of Globular Clusters display some of the most challenging cases 
of crowding that current instruments are capable of dealing with. 

The preparation of the simulated data can be summarized as 
follows: We used the V- and / band photometry of the glob- 
ular c l uster 4 7Tuc from|Sarajedini et al. (2007]) and |Anderson| 



|et al.| ( [20 08). To assign a realistic spectrum to each star, we 
constructed a s ingle isochrone (t = 13 Gyr, Z = 0.0045), us- 



ing the tool of Marigo et al. (2008] 



For each star we thus 
obtained estimates of and \ogg that were in turn used to 
select an appropriate spectrum from the stellar library of Munari 
et al.| ( [2005] ). The library spectra have an intrinsic spectral reso- 
lution of R = 20000 and cover the entire wavelength range from 
2500 A to 10000 A. In order to resemble our observational data, 
we extracted only the region around the Calcium triplet and con- 
volved the spectra with a Gaussian to obtain a final resolution of 
R - 7000. 

To simulate realistic datacubes, random fields of 8" x 8" (the 
FoV of the PMAS instrument) were selected from the central 
region of 47Tuc. Stars in those fields were represented by their 
respective spectra, and a velocity drawn randomly from a normal 
distribution with cr = lOkm/s, similar to the velocity dispersion 
in the centre of this cluster ( [McLaughlin et al. 2006 ), was as- 
signed to each spectrum. Using an analytical PSF profile with a 
full width at half maximum (FWHM) corresponding to a seeing 
around 1.0", cubes with 16 x 16 spatial pixels were prepared. 
Finally, appropriate noise was added. The resulting spatial sam- 
pling was 0.5"/pixel, again to resemble observational data ob- 
tained with PMAS. An example of the simulated data is given 
in Fig. [2] In the following, we refer to those datacubes as the 
'crowded field datacubes'. 

Besides these rather realistic datacubes, we also prepared 
idealized cubes containing only two stars. We used these to in- 
vestigate different aspects of the deblending procedure in iso- 
lation. The spectra of the stars were assigned using the same 
procedure as for the crowded field datacubes. We refer to these 
simulated data as the 'two star datacubes'. 



4. Deblending and extraction of spectra: Algorithm 

4.1. Global model 

We first clarify our adopted symbol and naming scheme. The 
datacube at hand is supposed to have pixel values b ;j ^, where 
the first two indices relate to spatial coordinates and the index 
k to the spectral axis. To denote different stellar sources in the 
datacube we use the superscript n while background components 
have a superscript m. A model datacube is then described as the 
sum of all sources that contribute flux, 



(1) 



In Eq. [T] f£ is the total monochromatic flux of star n in layer 
k, and psffj k is the fraction of f£ in the considered pixel. This 
fraction is equal to the value of the normalized PSF of star n 
at the spatial position (/, j) and at spectral pixel k. The flux in 
each pixel of background component m is given by bPj k . For a 
spatially constant background (such as the night sky), b™j k will 
only depend on k. An example for a background whose intensity 
varies across the field of view is the contribution of unresolved 



available at http : //stev . oapd . inaf.it/cgi-bin/cmd_2 . 3 
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Fig. 1. Impressions from the observations of several globular clusters that we obtained using the integral field spectrographs PMAS and ARGUS. 
The upper row shows a 20" X 20" cut-out from an HST/ACS image with the footprint of a single integral field observation. A whitelight image of 
each integral field data cube is shown below. In each image, north is up and east is left. 



stars to the observed flux; this will be discussed in Sect |4.6| The 
second sum in Eq.[T]can be further expanded to also include spa- 
tially resolved sources that might be present in the datacube. An 
example for such a component might be gaseous line emission 
in the observation of a young star cluster. Here we restrict the 
discussion to fields in where the flux is dominated by stellar con- 
tinuum sources. 

Starting from Eq. [T] we then search for the best (in a least- 
squares sense) model for a given dataset. We therefore have to 
minimize 



X 



E 

hj,k 



(2) 



Here, crf. k is the variance tailored to each pixel value. Un- 
fortunately, straightforward minimization of Eq. [2] is computa- 
tionally very demanding for the large parameter space that needs 
to be covered: in each layer k, a stellar source will contribute 
3 free parameters (/^ , x n k and z/jj) and a background component 
will contribute at least one additional free parameter. In addi- 
tion, several free parameters might be needed to find a suitable 
model for the PSF of the observation. A further complication is 
that Eq. [^represents a non-linear minimization problem. 

In order to make the search for a solution feasible, we split 
the optimization into three tasks: (i) an optimization for the PSF, 



(ii) an optimization for the source coordinates, and (iii) an op- 
timization for the fluxes. In each step, the model is optimized 
only with respect to one of these three properties, while the other 
two remain fixed to their current value. After one step has con- 
verged, the model parameters currently in focus are updated to 
their best-fit values and the model is optimized for the next set 
of parameters. On each layer k, the steps (i) to (iii) are then iter- 
atively repeated until convergence is found for the fluxes of the 
sources. 

The practical implementation of our approach can be sum- 
marized as follows: we start with an initial guess for the PSF and 
the source coordinates and use them to fit the fluxes in the cen- 
tral layer of the datacube. The reason for starting at the central 
layer is that the efficiency of the spectrograph should be high- 
est at the center of the covered wavelength range and therefore 
the data should give the tightest constraints on the model. After 
the iteration has converged to a solution for the central layer, we 
continue with the two adjacent layers. In this way, the analysis 
proceeds simultaneously to the red and blue end of the covered 
wavelength range. An integral field datacube provides a con- 
venient structure for such an iterative approach, as the changes 
in the PSF or source coordinates between two adjacent layers 
are always small. So the best-fit model of the previous layer is 
already a very good starting point for the analysis of the next 
one. A potential drawback of this approach is that if a single 
layer returns a strongly deviating model, such caused by an un- 
detected cosmic-ray hit or strong telluric absorption, it will affect 
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Fig. 2. Illustration of how the simulated data was produced and processed. Left panel: Input colour-magnitude diagram of 47Tuc. Overplotted 
is the isochrone that was used to assign spectra to the individual stars. Upper panels: Cut-out from an HST image of 47Tuc (upper middle); 
whitelight image of a simulated crowded field datacube centred on the HST footprint (upper right). The lower right panel shows the deblended 
spectrum of the star marked with a white circle. The input spectrum of this star is overplotted. 



the analyses of following layers. To avoid this, we generate ini- 
tial guesses by averaging over the best-fit models of the last N 
layers, with N typically being of the order of 10. 

In general, we do not consider cosmic ray hits as a major 
problem. The structure of raw IFS data actually allows for a 
convenient method to remove cosmic rays before performing the 
data reduction (Husema nn et al.|2 012). Additionally, it is pos- 
sible to significantly reduce the influence of undetected cosmic 
rays in the analysis. |Stetson| (p"987) presented a robust scheme 
that dynamically reduces the weights of pixels with strong resid- 
uals that is also applicable to IFS data. 

After all layers of the cube have been processed, their re- 
spective best-fit models are combined to produce a coherent 
wavelength-dependent PSF model and to determine source co- 
ordinates as a smooth function of wavelength. Using this infor- 
mation, the datacube is then processed for a second time, yet this 
time only the fluxes are fitted and we obtain the final spectra of 
the sources. 

Note that after a few layers have been analysed, the analyses 
of the red and blue half of the datacube proceed independently 
from each other. Thus, comparing the results obtained at the red 
and blue end can be used to check the reliability of the obtained 
model. 



4.2. Modelling the PSF 

The usual approach to determine the PSF in crowded field pho- 
tometry is to select a number of relatively isolated stars and fit 
them with an analytical function. To account for possible mis- 
matches between the analytical profile and the shape of the true 



PSF, an empirical look-up table correction is frequently applied 
afterwards. 

As discussed above, this approach cannot be applied to IFS 
datacubes without modification, mainly because of the very 
small FoV. Our adopted approach is to instead use the full en- 
semble of resolved stars within a field to reconstruct a global 
PSF model, obtained by means of a least squares fit of Eq. [2] 
to the data. The approach to recover the PSF using all stars in 
the field has previously been used by Schech ter et al.| ( [T993] ). A 
notable difference of our implementation is that the PSF is fit- 
ted to all stars simultaneously instead of in a sequential manner. 
Similar to Schechte r et al.| ([1993 ), we restrict the PSF descrip- 
tion to a purely analytical model, since the construction of a re- 
liable look-up table not only requires sufficiently isolated stars, 
but also a very well- sampled PSF. With the somewhat coarse 
spatial sampling of many IFUs this cannot always be taken for 
granted. But as the typical signal-to-noise ratio (S/N) per spec- 
tral pixel is much lower than in broad-band images, an analytical 
model should be adequate in most cases. 

To define an analytical PSF model we follow the approach 
adopted in GALFIT ( |Peng et al.(2002| ). In general, the PSF will 
not be round but rather of elliptical shape, with ellipticity e and 
position angle 0. To account for the ellipticity, the pixel coordi- 
nates x and y are first transformed into a coordinate system (x, y) 
centered on the origin of the star and whose x-axis is aligned with 
the semi-major axis of the PSF: 



x = (x - x n ) cos 6 - (y - y n ) sin #, 



y = ( x - x n ) sin + (y- y n ) cos 6, 



(3) 



(4) 
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Fig. 3. Radial surface brightness profile of the star BD+33D2642 
as measured in a PMAS datacube. Small grey dots show the flux in 
individual pixels, the black squares were obtained by radially binning 
the data. We used a Moffat profile and a simple Gaussian profile as the 
analytical PSF description and fitted the star with both PSF models. The 
blue solid line shows the best-fit PSF when using the Moffat, the green 
dashed line shows the best-fit PSF for the Gaussian profile. 



with (x n ,y n ) being the pixel coordinates of source n. The dis- 
tance to the centre of the PSF can now be written as as an angle- 
dependent quantity, 



r(x, y) = 




(5) 



In order to describe the radial shape of the PSF we adopt the 
Moffat profile, with a functional form given as 



M(x 9 y) = Jk 



1 + 



r(x, y) 



(6) 



The width of the Moffat profile is mainly determined by the dis- 
persion radius r d , while the /^-parameter defines the kurtosis of 
the profile, i.e the broadness of the wings of the PSF. The FWHM 
of the Moffat profile can be expressed in terms of r d and ft as 



FWHM = 2V2W- lr d . 



(7) 



This leaves us with 4 free PSF shape parameters per layer: ft, 
FWHM, e, and 6. Depending on the quality of the data to be 
analyzed, the number of free parameters can be reduced, for ex- 
ample by assuming a Gaussian instead of a Moffat profile, or by 
enforcing a circular PSF. The central intensity So of the PSF is 
directly tied to the monochromatic flux of a source. 

For each star n in the analysis, an empty array of the same 
size as a layer of the datacube is created. Its pixel values are set 
to the intensity of a normalized PSF at a radius r(x, y), using the 
current best-fit values for the source coordinates {x n v y n k ). Yet di- 
rectly using Eq.[6]to determine the pixel values psf"^ can lead to 
systematic errors close to the origin of the profile if the sampling 
of the PSF approaches critical values. In such cases, the change 
of M(x, y) within one pixel is so strong that its value at the cen- 
tre of a pixel is not a good approximation for the integrated PSF 
intensity in that pixel. Currently, we solve this problem by super- 
sampling each pixel within a certain distance to the centre of the 
PSF by a factor of typically 25-100 and calculate the intensity 
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Fig. 4. Example for the recovery of the PSF in a datacube. Shown 
is the accuracy of the recovered PSF parameters p {top) and FWHM 
{bottom) in comparison to their true values (blue dashed line). The grey 
solid line gives the best-fit values of the PSF parameters obtained in 
each layer and the final parameters of the PSF model obtained from 
these values are shown as a solid green line. Note that the analysis was 
started at the central image of the datacube and proceeded simultane- 
ously to the red and blue end of the cube. 



for each subpixel. The final value of each pixel then is obtained 
by summing over its subpixels. Numerical integration schemes 
such as the one presented by Buonanno & Iann icola| ( |1989| ) can 
significantly reduce the number of required subpixels and will be 
considered in our ongoing development of the algorithm. Once 
a PSF has been prepared for each source, we can use Eq|2] to 
optimize the PSF model for each layer of the datacube. After all 
layers have been processed, the results of the individual layers 
are modelled as a smooth function of wavelength for each free 
parameter of the PSF model. Following Wisotzki et al. ( 2003), 
we use low order polynomials for this task. This way we finally 
obtain the wavelength dependent PSF. 

To illustrate that our approach results in a valid description 
of the PSF in integral field data we analyzed one of our PMAS 
datacubes of the standard star BD+33D2642. The datacube con- 
tains a single star and thus allows for a precise measurement 
of the PSF even in the faint wings. We analysed the data in 
the way just described, using a single point source and a back- 
ground component. In Fig. [5] a comparison between a PSF as it 
is measured from our PMAS data and our analytical profiles is 
shown. It is obvious that a Moffat profile provides a good overall 
representation of the PSF. On the other hand, using a Gaussian 
severely underestimates the wings of the PSF. 

We used the simulated 'crowded field datacubes' to investi- 
gate how well our approach can recover the PSF in a dense stellar 
field. We produced 100 cubes based on a wavelength-dependent 
PSF represented by a Moffat profile with constant J3 and varying 
FWHM. An example for the PSF recovery is shown in Fig. [4] 
It demonstrates that the constraints on the PSF in an individual 
layer are not very stringent, as can be seen by the large scatter 
especially for the J3 fits. Yet the final combination of all individ- 
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Fig. 5. Comparison between recovered and input PSF parameters for 100 simulated crowded field datacubes. The PSF was modelled using 
a Moffat profile with wavelength-dependent FWHM and constant J3 parameter. Shown is the mean fractional difference in percent, between 
recovered and true FWHM (left) and p parameter (right), respectively, as a function of the input values. A single outlier of A/3 ~ 60% falls out of 
the plotted range and is shown as a lower limit. 



ual fits into a wavelength-dependent model yields a very close 
representation of the true PSF. 

To obtain a statistical impression on the accuracy of the PSF 
recovery we analysed the full sample of 100 simulated 'crowded 
field datacubes' and compared the recovered PSF parameters to 
the input. The results of this comparison are shown in Fig. [5] 

These plots show that the FWHM of the PSF can usually be 
recovered to an accuracy of 1-2%, while the recovered value of 
P usually agrees with the input value to within 10%. Thus, the 
PSF width is better constrained than the J3 parameter, similar to 
what we observed in Fig. [4] This is well known and due to the 
fact that /3 governs the outer wings of the PSF where the signal 
is generally low. But this also implies that errors in J3 do not 
influence the overall optimisation (Eq.[2]) as much as mismatches 
in FWHM. 

Furthermore, it is reassuring that the recovered values of the 
FWHM are unbiased in the sense that they scatter around the 
'true' value. This is not strictly the case for (3, which tends to 
come out slightly too high (by a few percent), thus correspond- 
ing to a PSF with slightly less pronounced wings than was put 
into the simulations. This is probably due to the fact that for 
high source densities, the ensemble of faint wings of each PSF 
can be misidentified as being part of the background; so appar- 
ently a very small fraction is, on average, transferred out of the 
wings of the stars used for the PSF estimation and into the back- 
ground component. We discuss the consequences of the achiev- 
able PSF accuracy on the quality of the extracted spectra in detail 
in Sect.0 

4.3. Source positions 

We assume that a source catalogue based on high-resolution 
imaging already exists. The task is thus to find a global trans- 
formation from a reference system (i.e. the source coordinates 
in the catalogue) to the coordinate system of the datacube. At 
least four parameters are required to describe this transforma- 
tion: a rotation angle a, a pixel scaling factor £ and a shift along 
both spatial axes, C and D. If we denote the coordinates of an 
object in the reference system by (u n , v n ), we can write down the 



coordinate transformation for a single layer as 

x\ = (cos orjfc u n + sin a k v n ) + , 
y\ = (- sin a k u n + cos a k v n ) + D k . 



(8) 
(9) 



Substituting A k = & cos or^ and B k = & sino^, Eqs. [8] and [9] can 
be rewritten as 



x n k =A k u n + B k v n + C k , 
y n k =A k v n -B k u n +D k . 



(10) 
(11) 



Eqs. 10 and 11 define a system of 2N linear equations, where N 
is the number of sources taken into account. While the transfor- 
mation will be wavelength dependent, this dependency can be 
further constrained taking into account that: 

(i) for instruments such as PMAS or ARGUS, where each 
spaxel is coupled to an optical fibre guiding the light into the 
spectrograph, neither a nor £ should depend on wavelength and 

(ii) the variation of both C and D with wavelength due to 
atmospheric dispersion can be predicted using a model for the 
wavelength dependent atmospheric refractive index during the 
observation. Such predictions as a function of airmass and par- 
all actic angle have been derived, e.g., by Filippenko (1982]); also 
see Sandin et al. ( 2012 ) for an overview of different approaches. 
For the range of airmasses (< 2) of our PMAS data, an offset 
up to 072 (equivalent to 0.4 spaxels) is predicted across the cov- 
ered wavelength range, quite similar to the values we actually 
measure. 

One might use this information to eliminate the wavelength 
dependency in Eqs. [TO] and [TT] already in the data reduction and 
then assume a single transformation for all layers of the cube. 
However, such a correction would involve resampling the data 
which we believe is to be avoided as much as possible. Instead, 
we start by finding a best-fit solution for every layer separately 
and then use the information to constrain the transformation a 
posteriori: After all layers have been analysed, we determine the 
IFU coordinates x n and y n of every source as a smooth function 
of wavelength. To fix the rotation angle and the pixel scaling 
to a common value for all layers of the datacube, the polyno- 
mial fits can be coupled in a way that the change in x n and y n 
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Fig. 6. Comparison between the recovered and true source coordinates for the sources in 100 simulated crowded field datacubes, expressed as 
the mean deviation in x- and ^-direction relative to the FWHM of the PSF as a function of the position of each source, in units of the FoV size. 
Dashed vertical lines indicate the FoV edges. The median deviation along x and y is marked by a solid black line, the 68.3% and 95.4% percentiles 
of both distributions are indicated by the dotted and dashed black lines, respectively. A few outliers are out of the plotted range and indicated by 
upper limits. 



with wavelength is the same for all sources n. Furthermore, that 
slope can be fixed to the one predicted by atmospheric disper- 
sion. This approach is very flexible with respect to how much a 
priori information should be used in determining the coordinates 
(e.g., the determination of the refractive index depends on an ac- 
curate knowledge of atmospheric properties such as humidity or 
air pressure that is not always available) and it uses the informa- 
tion from all layers together to obtain a final result. 



Note that instead of fitting the parameters of Eqs. 10 and 
[TT] we directly fit a low order polynomial to each x n and y n . 
This is done because the variation of A to D with wavelength 
is not completely independent from one another. For example, 
to a certain extent a change in the scale factor £ (and thus A 
and B) might be compensated by changes in the shifts C and 
D, especially when the S/N in a layer is low. Such correlated 
variations of the parameters on small wavelength scales would 
be washed out by a polynomial fit, causing larger errors in the 
recovered IFU coordinates. 



In this simple linear model, Eqs. 10 and 11 do not account 
for distortions of the FoV. Given the small number of spatial 
pixels provided by most IFUs, this approximation should always 
be valid. If that should be not the case, a solution would be 
to expand the coordinate transformation to also include higher 
order terms of u n and if. This will likely become important for 
more complex instruments such as MUSE (cf. Sect[7]), where the 
FoV is split into 24 different parts that follow independent opti- 
cal paths through the instrument before finally getting dispersed 
by different spectrographs. 

One important feature in our procedure is that the global 
transformation model can include sources with centroids actu- 
ally outside of the FoV. Of course the accuracy of the positions 
will decrease with increasing distance since the transformation 
itself can only be constrained inside the FoV, but we are mainly 
interested in sources close to the observed field which still can 
have an influence on the light distribution inside the FoV, espe- 
cially in the case of bright stars. 

Again, we used our 100 simulated crowded field datacubes 
to assess the accuracy of the recovered positions. In Fig. [6] we 
plot the offsets between the input and the recovered source posi- 



tions, in units of the PSF FWHM in each cube, so that the results 
provide a generic measure of the achieved accuracy relative to 
the achieved spatial resolution. 

We find that the standard deviations of the recovered coordi- 
nates from the true ones are of the order of 1-2% of the FWHM 
of the PSF. The highest accuracy is achieved in the centre of 
the FoV where the coordinate transformation is best constrained. 
The scatter increases only slightly towards the edges of the FoV, 
but there is also a small systematic offset in the sense that at 
small x- or ^-coordinates, the recovered values are on average 
smaller than the true ones, while they are larger on average for 
large coordinates. With regard to the coordinate transformation, 
this implies that the scale factor of the transformation & is re- 
covered too large. This behaviour can be attributed to the small 
FoV that we have to deal with. It leads to a significant fraction 
of sources that contribute to the observed flux distribution hav- 
ing their centres outside the FoV. In case of a small mismatch 
between the recovered and the true position, the strongest resid- 
uals will emerge in the region around those two positions. The 
more this region is pushed away from the observed field, the 
smaller its effect on x 2 will t> e - Thus x 2 (£) is asymmetric to- 
wards higher values of causing the observed trend. One can 
obviously avoid this behaviour by making assumptions about the 
coordinate transformation and fixing the value of However, in 
this case the pixel scale of the integral field spectrograph must be 
known very precisely, to better than 0.01 " to achieve a compa- 
rable accuracy. We emphasize that this trend is extremely small 
and that in many cases the uncertainties of the individual high 
resolution coordinates will already be higher. 

4.4. Extraction of spectra 

Once the PSF model and all source positions are known with suf- 
ficient accuracy, the least- square solution for the source spectra 
becomes a linear equation. We have to minimise 

* 2 = |Aa-b| 2 . (12) 

We denote A the 'PSF matrix' because it contains the PSF of 
every source in the fit. a contains the object fluxes we aim 
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at and b is the data. Note that, because of the 3 -dimensional 
data structure, a and b are also 3 -dimensional and A is even 4- 
dimensional. We can then express the individual components of 
Eq.[T2|as 



A n - 

A i,j,k ~ 



^Uk 
°~ij,k 

r, 

hj,k 



i,j,k 



(13) 
(14) 
(15) 



If we want to solve Eq. [12] in a least- squares sense we have 
to properly take into account the uncertainty cr,-^ of each data 
value. To achieve this, both the data and the PSF coefficients 
have to be divided by the uncertainties ( [Press et al.|1992| ). It can 



15 into Eq. 12 



easily be verified that substitution of Eqs. [13) to 
yields Eq.|2] 

The solution to Eq. [12] is obtained separately in each layer 
of the datacube without any coupling of the fluxes of an in- 
dividual star in adjacent layers. In this sense, it is similar to 
performing photometry on each star in each layer (with known 
PSF and source coordinates) and obtaining the spectrum of each 
star as the combination of its individual monochromatic fluxes. 
However, in contrast to common crowded field photometry, the 
monochromatic fluxes in one layer are obtained simultaneously 
for all stars instead of measuring star after star. 

For each spectrum that is obtained by solving Eq.[l2]we can 
estimate the S/N per layer via 



S/N k = f k 



( 

E 



P Sf ^ 



xl/2 



'Uk 



(16) 



Note that Eq. [16] only takes into account the noise expected due 
to the pixel uncertainties cr l and therefore gives only an upper 
limit which does not include any noise contributions introduced 
by connecting the pixels in the analysis. As we will show below, 
the true S/N in a deblended spectrum can be significantly lower. 

It is worth spending a few more words on the handling of un- 
certainties in the analysis. Throughout this paper we assume that 
the uncertainties are known and correct, i.e that the true variance 
in each pixel is cr?^. The reduction of integral field data usually 
requires at least one step of resampling the data onto a regular 
grid. Resampling always introduces covariances between neigh- 
bouring pixels, although they are usually neglected in the further 
analysis. Consequently, the provided variance tends to under- 
estimate the true uncertainty of a pixel value. However, this has 
no effect on the optimization which is based on^ 2 minimization. 
Only the interpretation of absolute^ 2 values with respect to con- 
sistency between a model and the data may become doubtful. 

Since Eq.[T2] describes a linear least squares problem, its so- 
lution does not require any starting values and can be directly 
obtained by matrix inversion. For a realistic number of a sources 
this is not a computationally intensive process, even if applied to 
all layers in a cube. Furthermore, the computing time required 
for the solution of Eq.[l2]can be very significantly decreased fur- 
ther by taking into account the fact that the contribution of any 
star is limited to the pixels in its vicinity. Each element n of the 
PSF matrix A therefore contains a large number of pixels with 
essentially zero values. For this reason, A can be written as a 
sparse matrix, for which there are dedicated efficient algorithms 
available (e.g., Paige & Saunders 1982). With these provisions 
it becomes possible to fit the fluxes of even several thousand 



sources in a MUSE datacube simultaneously (see Sect.[7]below). 
The computation time required to analyse a datacube is almost 
exclusively determined by the time required to obtain the PSF 
and the coordinate transformation. The computation time scales 
linearly with the number of sources. When fitting ~ 20 sources 
for an instrument comparable to PMAS, one iteration on a single 
layer takes around 10 s on a single CPU. For the typical number 
of iterations and layers in a cube, the total time required for a 
whole datacube is around 10-20 hours on a single CPU. Some 
parts of the code have already been parallelized, others will be 
in the future, so that the actual time required for the analysis can 
be strongly reduced. 

4.5. Construction of the source list 

When the input catalogue of sources in the field is constructed 
from high-resolution imaging, it will typically reach magnitudes 
where the S/N level in the IFS data is too low to produce mean- 
ingful spectra. We thus need to construct a subsample of stars 
whose spectra can possibly be deblended in an available dat- 
acube. By design, this subsample will contain stars over a large 
range in magnitudes and expected S/N levels. Even if we are 
interested mainly in the stars bright enough to yield a spectrum 
with a sufficient S/N for some analysis, we also need to account 
for the effects of blending with the (more numerous) fainter ob- 
jects. There is a limit, however: When true source confusion 
sets in, the deblending process and the flux assignment to indi- 
vidual sources becomes to some extent arbitrary. In Sect. [5] we 
explore quantitatively by means of simulations where this limit 
is reached. 

The decision of whether or not to include a particular source 
will depend on several criteria: (i) the brightness of the source, 
(ii) the distance to other nearby sources and their relative bright- 
nesses, (iii) the position of the source, in particular if it is located 
close to the edge of the FoV (or even outside, see below). Effec- 
tively, the first two criteria can be combined into a single one 
based on S/N, taking the degradation of S/N due to crowding 
into account. 

The practical sequence of constructing a source list is as fol- 
lows. We first estimate a limiting magnitude where confusion 
becomes dominant, based on a global characterisation of the ex- 
posure, given its depth, resolution etc., in comparison with the 
input catalogue. We then select a preliminary source list on the 
condition that the source magnitudes are brighter than the con- 
fusion limit. For those sources we predict the continuum S/N 
using simulations as explained below, accounting for the overall 
effects of crowding as well as for the influence of close-by bright 
stars. The final source list is then based only on the expected S/N 
of the spectra. 

In this process there will be stars with magnitudes brighter 
than the confusion limit that do not pass the second selection 
stage, given their proximity to a brighter star. Yet the influence 
of those stars has to be taken into account in the analysis. In such 
cases we generate a modified PSF for the close-by bright star 
that approximates the contribution of its fainter companion using 
the broad-band magnitudes of the two stars and their distance. 
The extracted spectrum will then be a combined one, and will be 
flagged accordingly in the resulting catalogue of spectra. 

We note that the selection using S/N ratios will also take care 
of picking the sources close to the edges of the FoV that con- 
tribute significantly to the observed data. Eq. [16] is a sum over 
all spatial pixels, the flux that enters in the S/N calculation is the 
fraction of flux that is recorded by the detector. For a given posi- 
tion, this fraction is determined by the PSF. The further a source 
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is located outside the FoV, the smaller this fraction, and thus the 
S/N will be. 



4.6. Treatment of background 

So far we concentrated our discussion on the resolved sources. 
But the data will always contain some flux contribution from 
unresolved background components. One such contribution is 
the night sky, the intensity of which we assume to be constant 
over the field of view. Just like any other flux component that 
is spatially flat, it can be easily accounted for by expanding the 
PSF matrix A by a component whose values are inversely pro- 
portional to the uncertainties cr^^ 

Of greater interest is another background component that is 
produced by stars below the confusion limit, i. e. th ose that are 
not included in our source list discussed in Sect. 14.51 This back- 
ground will in general not be spatially flat, but instead the distri- 
bution of the stars will produce a grainy structure, known from 
the surface brightness fluctuations observed in nearby galax- 
ies. Our approach allows us to actually model these fluctua- 
tions. The sources above the confusion limit provide us with 
a model of the PSF in the observation. Additionally, from the 
high-resolution imaging in combination with our model for the 
coordinate transformation discussed in Sect. |4.3| we obtain pre- 
cise IFU coordinates even for the sources below the confusion 
limit. Furthermore, we can use the photometry provided by the 
high-resolution imaging to estimate the relative brightness of the 
sources. With PSF, positions and relative brightnesses in hand, 
we can predict the relative brightness of the grainy background 
in every spatial pixel of the datacube. This prediction can then 
be included (again after division by cr ZJ ^) in the PSF matrix. 

Another physical background component could be in the 
form of gaseous emission, e.g., from filaments in an Hn region. 
Obviously, such a component is not to be expected in globular 
clusters but we want to provide a tool that is useful in any kind 
of crowded stellar field. We also mention again that the scien- 
tific focus might actually not lie on the resolved stars but the 
nebular emission: One might be interested in the kinematics or 
emission line ratios and wish to remove the contamination of 
the bright stars. In their work on extragalactic planetary nebu- 
lae, jRothetaL] ([2004]) showed that integral field spectroscopy is 
capable of separating individual point sources from diffuse emis- 
sion of the ISM. The ISM component will be line emission and 
thus be restricted to a few layers in the datacube. So even if it 
strongly biases the PSF fit in those layers, a reliable PSF model 
can be obtained by interpolation of the results obtained blue- 
wards and redwards of the emission line. One might then start 
from some initial guess for the spatial intensity distribution of 
the gaseous emission that is included in Eq. [T2] and iteratively 
improve it based on the residuals observed in the layers under 
consideration. 



5. Performance of the deblending process 

The quality of the spectrum extracted from a datacube contain- 
ing a single isolated star will depend almost exclusively on the 
brightness of that star. This is different for the extraction from 
crowded fields, where several effects may contribute to degrade 
the spectrum. We investigated the performance of our deblend- 
ing and extraction algorithm on the basis of our simulated data- 
cubes. We employed several criteria to quantitatively measure 
the quality of the extracted spectra: 
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Fig. 7. Degradation of the S/N ratio in the spectrum of a star with 
a nearby neighbour, as a function of the distance to the neighbouring 
source (normalized to the width of the PSF) and of the magnitude dif- 
ference (star minus neighbour) between the two stars. We define the 
relative S/N as the ratio between the measured value and the value ex- 
pected for an isolated star. Note that for clarity, small horizontal offsets 
have been applied to the plotted datapoints for the different magnitude 
differences. 



The first two criteria are purely formal indicators: (i) We 
measured how the S/N in the continuum behaves under crowd- 
ing conditions, relative to the value in isolated stars of the same 
magnitude, (ii) We also tested how robustly the continuum level 
is recovered, in terms of a potential systematic error in the broad 
band magnitude. 

In order to facilitate a discussion of the astrophysical pos- 
sibilities and limitations of crowded field spectroscopy, we also 
considered the behaviour of two types of derived spectral pa- 
rameters on crowding: (iii) radial velocities, and (iv) equivalent 
widths of strong absorption lines. 

In the following we illustrate and discuss the performance 
of our code for each of these parameters. We first illustrate and 
quantify crowding effects for the restricted scenario of only two 
stars within a cube, with varying angular distance and brightness 
difference. We then explore the global performance in realistic 
situations using the 'crowded field datacubes' representing mock 
PMAS observations of a real globular cluster, constructed as de- 
scribed in Sect. 13.21 



5.1. "Crowding" of two stars 

5.1 .1 . Signal to noise ratio and continuum level 

To assess how the S/N is affected by crowding, we compared 
the S/N expected using Eq. [16] to the one actually measured in 
the deblended spectra. To measure the S/N of a deblended spec- 
trum we did the following: we subtracted the (noise free) input 
spectrum from the deblended one. To account for possible mis- 
matches in the continuum level and slope of the extracted spec- 
trum we fitted the residuals with a low order polynomial that was 
then also subtracted. The residuals should now scatter around 
zero with a standard deviation equal to the noise level of the 
deblended spectrum. We determined the standard deviation in 
a window around the central wavelength and divided the mean 
flux of the deblended spectrum in that window by it to yield a 
S/N ratio. 

For two stars of varying angular separation and magnitude 
difference, Fig. [7] shows the dependence of the S/N degradation 
on source separation d n and flux ratio. For separations greater 
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Fig. 8. The average error in the continuum level of the deblended 
spectra (expressed as a magnitude difference) as a function of the dis- 
tance of the source to a neighbour (again normalized to the width of the 
PSF) and its relative brightness. As in Fig. [7] datapoints for different 
flux ratios are plotted with small horizontal offsets for clarity. 



than the FWHM of the PSF, no degradation occurs as expected. 
For smaller separations the S/N decreases steadily and is only 
half that of an isolated star for dn - 0.3 FWHM. The reason 

" As 



for this degradation is an increasing degeneracy in Eq. 13 



the locations of the two stars approach each other, their PSF im- 
ages become nearly the same, and the solution of Eq. [13] is no 



longer unique; any linear combination of the two stars that main- 
tains the total flux will yield an almost equally likely (in terms of 
X 2 ) result. Practically this means that in every layer a different 
amount of flux will be transferred from one star to the other. 

An interesting aspect revealed by Fig. [7] is that the decline 
in S/N seems to be almost independent of the flux ratio between 
the two stars. This is unexpected at first sight - one might expect 
a bright companion to have a higher impact on the deblended 
spectrum of a faint star than vice versa. But in fact the observed 
behaviour follows directly from Eq. [16} For small separations, 



the expected noise (the inverse of the square root in Eq. [16]) will 
be almost the same for the two sources. But the final noise will 
also be the same for both sources because it is determined by the 
amount of flux that is shuffled back and forth between the two. 
Of course, the absolute S/N will still be higher for the brighter 
source because its signal is higher. 

It is interesting to view the results depicted in Fig. [7] in 
the light of the Rayleigh criterion, which states that two point 
sources are resolved if their separation is larger than the FWHM 
of the PSF. It is of course well known that this criterion is not 
a strict limit but only indicative for the spatial resolving power. 
Fig. [7] shows that in the case where the relative positions of the 
two sources are known, it is possible to deblend stars with dis- 
tances well below the FWHM of the PSF, although at the price 
of a reduced S/N. 

We also investigated how well the continuum level of the 
stars could be recovered. To this aim we characterise the contin- 
uum by dividing each deblended spectrum by the corresponding 
input spectrum and converting the result into a magnitude. In 
Fig. [8] the dependence of the average continuum error on source 
separation and flux ratio is shown. It is remarkable that the 
strong decrease in S/N visible in Fig. [7] does not cause a simi- 
lar degradation of the continuum level. This implies that the in- 
crease in noise is indeed purely random and does not introduce 
any systematics in the deblended spectra. The fact that the actual 
deblending of the spectra (after PSF and source positions have 
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Fig. 9. Accuracy in the recovered radial velocities {bottom) and 
equivalent widths of strong absorption lines (top) for a blended star, as 
a function of source separation and for different flux ratios. The lines 
depict the median velocity difference (measured - true, thick solid) and 
the 75% percentiles (thin dashed) for a distribution of 1000 stellar pairs 
per separation value. 



been determined) does not require initial guesses makes it very 
robust against systematics. Especially in such cases discussed 
above where the two stars are very close to one another and the 
X 2 value becomes insensitive to the flux ratio between the two 
stars, the outcome of a fit that requires an initial guess would 
strongly depend on the value of that guess. 

We do observe a small overestimation of the continuum level 
in the case of a very bright neighbour. However, the systematic 
error stays below O.lmag and one should keep in mind that in 
this cases we are trying to measure the flux of a star that has a 
companion well inside the extent of the PSF that is at least 15x 
brighter. 

5.1 .2. Radial velocities and equivalent widths 

We now consider the recovery of astrophysical quantities. We 
determined radial velocities by cross-correlating the extracted 
spectra with the noise-free input spectra used to generate the 
datacube. The quantity of interest is then the 'measured minus 
true' velocity difference Ay. Recall that the spectral resolution 
adopted in the simulations was A/AA = 7000, thus correspond- 
ing to a velocity resolution of 42.9 km/s (FWHM). 

The outcome of running our deblending code on a large num- 
ber of simulated 'two star datacubes' (1000 realisations per sep- 
aration setup) is depicted in the lower panel of Fig. [9] The scat- 
ter of the velocity differences increases as the source separation 
decreases, which is of course expected as an immediate con- 
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sequence of the declining S/N. However, the median values of 
Av stay comfortingly close to zero, although the two stars have 
randomly assigned velocity differences of the order of 10 km/s 
(see Sect. |3.2| ). This is even so if the neighbouring star is much 
brighter than the analysed source. 

But Fig.[9]shows also clearly that the impact of partial blend- 
ing in terms of the scatter on Av is strongest if the nearby star is 
significantly brighter. At first sight this behaviour might be seen 
as different from what we observe for the S/N (cf. Fig. [7]), but it 
can be easily explained by the fact that the accuracy in the radial 
velocity measurement depends on the actual S/N, rather than on 
the S/N degradation. 

To check whether the deblending process might lead to a sys- 
tematic error in Av we divided our sample of stars into two sub- 
samples, one with all stars that were assigned a positive 'true' 
velocity relative with respect to their neighbours and one with 
all stars with negative 'true' velocity. For both subsamples we 
found that the recovered median Av values are statistically indis- 
tinguishable from zero. 

To investigate the behaviour of absorption line equivalent 
widths we focused on the Calcium triplet at /L18498, 8542, 
8662 A. The equivalent width (EWcaiO of this feature is widely 
used to estimate the metallicity of stars (e.g., Batta glia et al.| 
2008J, and it is therefore important to look into the integrity of 
this quantity under crowding conditions. 

The behaviour that we observe for the recovered values of 
EWcaT is very similar to that just described for the radial veloci- 
ties. The upper panel of Figure[9] shows the median deviation and 
75% percentiles between the equivalent widths measured in the 
deblended spectra and those measured in the input spectra. On 
average there is again no bias, but the scatter increases with de- 
creasing separation and with increasing brightness of the neigh- 
bouring star, as expected. 



5.2. Performance in realistic crowded fields 

In the last section, we showed how we can predict the expected 
S/N of a deblended spectrum under crowding (cf. Fig. [7]), us- 
ing the idealized case of only two stars. But when analysing a 
crowded stellar field, we have to take also another effect into ac- 
count. Below a certain magnitude, the confusion limit, the stellar 
density of similarly bright stars will be so high that they form a 
'pseudo background' . When this limit is reached, longer integra- 
tion times will not lead to an increase of the number of resolved 
sources, though the average S/N will still increase. 

In order to facilitate the following discussion, we first intro- 
duce the term resolution element as the area covered by a circle 
whose diameter is equal to the FWHM of the PSF. When deal- 
ing with source densities, it is quite useful to specify them as 
numbers of sources per resolution element because this measure 
is independent of the specific instrument characteristics (number 
of spaxels, spaxel size) and observing conditions (seeing). Note 
that when stating source densities in the following, we refer to 
the density of stars brighter than a given limit. 

Imaging studies are usually considered fairly complete down 
to source densities of 0. 1 stars per resolution element. Of course, 
there is no sharp cut between detected and undetected sources 
at this limit as some brighter sources will already remain unde- 
tected while some fainter ones will still be found. Our analysis 
is based on an existing inventory of sources, so there is no need 
for a source detection. Instead, we define a resolvable source as 
one that still improves the overall quality of the deblending pro- 
cess when it is included. Later, we will also discuss the subset 
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Fig. 10. The recovered continuum level of extracted spectra as a 
function of the number of stars in the source list, from the analysis of 
100 simulated crowded field datacubes. The number of stars is given per 
datacube in the bottom and per resolution element in the top label. A 
thick dashed line is used to indicate the median continuum error of the 
brightest 10 stars per cube, and dotted lines enclose the 75% percentiles 
of the distribution. 



of useful resolvable sources, which are those for which physical 
parameters can be recovered to a given accuracy. 

The aims we pursue in this section are twofold. First, we 
want to obtain a well-founded determination of the confusion 
limit, i.e., the transition from resolvable sources to unresolvable 
sources, and investigate the effects of selecting either too few or 
too many stars. This is an important aspect for the source selec- 
tion process that was presented in Sect. |4.5| Second, we want 
to verify whether the effects we discussed for the crowding of 
two stars (cf. Sect |5.1| ) can also be identified in realistic crowded 
stellar fields. 

We used the crowded field datacubes and tested how the 
number of stars included in the deblending process influences the 
results by including all stars in the deblending process brighter 
than a limiting magnitude m cut . m cut was varied over a range that 
corresponded to average source densities between 5 and >100 
stars per simulated datacube. In the following discussion, we 
will use two measures for the source density: besides the num- 
ber of sources per resolution element we also give the number of 
stars per datacube, since absolute numbers are quite intuitive. In 
our application cases, one crowded field datacube contains 256 
spatial pixels (reproducing the characteristics of the PMAS in- 
strument) and because the FWHM of the PSF is 2 pixels, the 
number of resolution elements per cube is -80. 

5.2.1. Continuum biases 

We first checked whether the continuum level of the stellar spec- 
tra becomes biased after extraction from a crowded stellar field. 
Again we converted the fraction between recovered and true 
spectrum into a magnitude. So in the case of systematic flux 
transfer to or from other stars due to source confusion, we ex- 
pect a non-zero offset Af. We performed the deblending experi- 
ment on 100 simulated datacubes and measured the distribution 
of Af for the set of the 10 brightest stars in each datacube, as 
a function of the total number of stars in the source list. Re- 
call that for each cube, the total number of stars in the source 
list was varied from a few to > 100 per field. The results are 
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Fig. 11. Estimated S/N of extracted spectra as a function of the num- 
ber of deblended stars per field or per resolution element, respectively, 
again using the 10 brightest stars per simulated cube. The thick solid 
line indicates the median of the distribution. As in Fig. [7] we show the 
ratio of the S/N measured in a deblended spectrum and that expected 
based on Eq.[T6| 



based on the 10 brightest stars per cube to allow for a fair com- 
parison between the individual simulations. Using all deblended 
stars would penalize the simulations in which many stars were 
included because stars with fainter magnitudes will on average 
have deblended a spectrum with a lower S/N. 

Figure [10] shows the median and 75% percentiles of Af and 
its dependence on the number of stars in the source list. The 
median value of Af is nearly always very close to zero, im- 
plying that on average, the continuum level remains essentially 
unbiased under these crowding conditions. Only a very minor 
systematic offset to fainter magnitudes of ~ 0.01 mag is ob- 
served. We expect that the missing flux is lost in the wings of the 
PSF and transferred into the background as it would agree with 
the fact that the PSF is recovered with slightly less pronounced 
wings (cf. Fig. [5]). 

It also becomes clear from Fig. [TO] that individual stars can 
show significant deviations in their continuum levels. The ac- 
curacy in the recovered continuum is actually worst when only 
the few brightest stars (< 10 stars per simulated crowded field 
datacube) are extracted because in that case stars falling just be- 
low the selection cut are still well resolved, yet their contribution 
is not accounted for during the deblending. m cut is significantly 
brighter than the confusion limit and thus the stars are extracted 
against a very inhomogeneous 'pseudo' -backgro und. This be- 
haviour is quite similar to that reported by Moehler & Sweigart 
2006) when performing multi-object spectroscopy of horizontal 



branch stars in the globular cluster NGC6388: the contribution 
of close-by stars could not be accounted for and, as suggested 
by the authors, did likely influence the results. Our results show 
that such problems are essentially avoided when applying PSF- 
fitting techniques on IFS datacubes. The recovery of the contin- 
uum level is significantly increased upon including more sources 
in the deblending process. 



5.2.2. S/N degradation 

The median relative S/N of the 10 brightest sources in each cube 
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Fig. 12. Accuracy of the equivalent widths (top) and radial velocities 
(bottom) determined from the recovered spectra of the brightest 10 stars 
per cube as a function of the number of sources included. Line types 
are as in Fig.fTO] 



trum decreases when the total number of sources taken into ac- 
count is increased. This decrease can only be due to source con- 
fusion; in Sect. |5.1| we discussed the consequences of deblend- 
ing sources with a small mutual distance. Now, if we increase 
the number of stars in the process without taking into account 
the environment of a selected source, it becomes more likely 
that stars are included whose mutual distance is close to or even 
smaller than the minimum distance at which the spatial resolu- 
tion of the data allows for a clean separation of the two. If one 
tries to deblend those sources nevertheless, this will lead to an in- 
crease in the noise of the extracted spectra. Yet the comparison 
with Fig. [TO] also shows that this flux reshuffling is essentially 
random, individual stars on average do not receive flux or loose 
flux. Only their spectra get somewhat noisier. 

In the source selection scheme that we have adopted, such a 
behaviour is avoided by including the expected S/N in the selec- 
tion process. For each source, the S/N is estimated using Eq.[T6 
and applying the correction found in Fig. [7] Then only sources 
above a threshold in S/N are considered resolvable whilst others 
are added to a brighter neighbour. 



is shown in Fig. 1 1 Apparently, the S/N in a given stellar spec- 



5.2.3. Recovery of spectral parameters 

To quantify the influence of the number of deblended sources 
on our ability to recover physical parameters from the spec- 
tra, Fig. [12] shows the accuracy of measured radial velocities 
and equivalent widths as a function of the number of deblended 
sources. 
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Fig. 13. The fraction of recovered sources is plotted as a function of 
the source density, using three different values for the required accuracy 
in the recovered radial velocity. The curves have been obtained using 
the same data that yielded the most accurate results in Fig. [12] 



When only the few brightest stars in each cube are de- 
blended, the true EWc a T are underestimated. This can be easily 
explained as a consequence of flux transfer from unaccounted 
fainter to brighter stars, because the EWc a T increases as one 
moves up the red giant branch. Once the source list accounts 
for those fainter stars, this bias disappears entirely. 

More interesting is the behaviour of the radial velocity accu- 
racy Ay, which has a broad but clearly defined minimum. This 
suggests that there is actually an optimal number of stars to be 
used in the deblending process. The reason for the occurrence of 
such a minimum is that two counteracting effects influence the 
results, as we have identified previously. Inclusion of too few 
stars leads to systematic errors because the fluxes of fainter (yet 
still resolved) stars are not accounted for. Inclusion of too many 
stars on the other hand causes a drop in the S/N of the deblended 
spectra that renders the determination of stellar parameters less 
accurate. 

The source density at which those two effects are best bal- 
anced is -30 stars per simulated cube, corresponding to 0.4 stars 
per resolution element. If we adopt this number as the confusion 
limit and compare to the value of 0. 1 stars per resolution element 
in crowded field photometry, we can quantify the improvement 
we get from using a pre-defined source catalogue instead of hav- 
ing to detect the sources in the datacube. Furthermore, we can 
use this value to give requirements on the quality of the input 
source catalogue in order to avoid being limited by the num- 
ber of sources it contains: the stellar density of detected sources 
should be at least 4x higher than the spatial resolution of the 
integral field data would allow for. Thus, the spatial resolution 
of the observation used to create the input catalogue must be at 
least twice the spatial resolution of the datacube. 

So far, we have counted all stars as resolvable that yielded 
on average better results when included in the deblending pro- 
cess. Clearly, not every single star of those will yield a useful 
spectrum. The subsample of useful resolvable stars will strongly 
depend on the science goals. To demonstrate this, we used the 
simulations that yielded the most accurate results, i.e., where 0.4 
stars per resolution element were deblended and determined the 
recovery fraction as a function of source density when request- 
ing different levels of accuracy in the recovered radial velocities. 



The recovery fractions we obtained are shown in Fig. 13 For 
the very stringent condition that uncertainties are <1 km/s, the 
completeness drops to 50% already at 0.1 sources per resolu- 
tion limit. On the other hand, under more relaxed conditions, 
the completeness drops to below 50% only at sources densities 
>0.2 sources per resolution limit. 

Finally, we note that the recovered velocities have a very 
small systematic offset of -0.2 km/s, independently of the num- 
ber of stars; for the present study this is however of no concern. 

5.3. Influence of crowding 

Our crowded field datacubes represent realistic integral field ob- 
servations of a Globular Cluster, with stellar densities typical for 
Globular Clusters. We now want to investigate the performance 
of our method in different regimes of crowding. Over a certain 
range in stellar density we expect a trade-off between the crowd- 
ing and the achievable depth: The more crowded the stellar field, 
the more our analysis will have to be restricted to the brightest 
stars. Yet with increasing stellar density, it will also get more 
challenging to deblend clean single object spectra at all because 
the contrast between the individual sources decreases. In the 
limiting case the stellar density will be so high already for the 
brightest sources that the stellar field is entirely unresolved. 

In this section we aim to quantify up to what amount of 
crowding our approach yields useful results before it breaks 
down. Quantifying this limit is important when making predic- 
tions about whether a stellar field can be accessed by means of 
crowded field 3D spectroscopy. This will not only be applica- 
ble to Globular Clusters but to any type of crowded stellar field. 
In nearby galaxies, for instance, projected stellar densities can 
significantly exceed those of a typical Globular Cluster. Thus it 
would be very helpful to know up to what source density we can 
still obtain good results. 

To test the influence of the crowding, we modified our sim- 
ulations of crowded field datacubes in the following way: we 
identified as a bright star every source in the catalogue with a vis- 
ible magnitude brighter than the Horizontal Branch (F606W < 
13.5 for 47Tuc), i.e we concentrated on the brightest giants. For 
each simulated datacube, we randomly picked stars from the cat- 
alogue and placed them in the datacube until a certain number of 
bright stars was reached. The number of bright stars picked var- 
ied between 4 and 400. The further processing of the cubes was 
then similar to the simulations described in Sect. 13.21 

In total, 200 datacubes were prepared that were all anal- 
ysed using our algorithm. To quantify its performance for a 
single cube we counted the number of bright stars whose de- 
blended spectra fulfilled an accuracy criterion. Two different 
accuracy criteria were used: an error in the recovered contin- 
uum of <0.1mag and an offset in the recovered radial velocity of 



<2kms _1 . In Fig. [141 we show the number of recovered sources 
as a function of the crowding. "Crowding" here is defined as the 
number of bright sources in a datacube. Furthermore, we again 
normalized the star counts by the number of resolution elements, 
for the reason mentioned above. 

both accuracy criteria yield compara- 



As Fig. 14 



shows, 

ble results: We observe that up to a crowding of 0.2 sources 
per resolution element, the number of accurately deblended 
sources increases approximately linear with the number of exist- 
ing sources. For a higher crowding of 0.2-1.0 sources per resolu- 
tion element, we observe a plateau with an average of ~ 0.2 accu- 
rately deblended sources per resolution element. If the crowding 
increases beyond 1 source per resolution element, the number of 
deblended sources that fulfil our accuracy criteria starts to de- 
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Fig. 14. The impact of crowding on the efficiency of the deblending algorithm. We show the number of useful deblended spectra as a function of 
the number of bright sources in the simulation. To obtain a measure of the crowding that does not depend on the specific simulation setup (seeing, 
size of the field of view) we have normalized the number of sources by the number of resolution elements. Two criteria were used to identify a 
usefully deblended spectrum: a magnitude offset < 0.1 mag {left) and an offset in radial velocity < 2km s _1 {right) between the input spectrum and 
the recovered one. Each plotted datapoint corresponds to the analysis of one mock datacube. In both panels, a grey solid line indicates the optimal 
case where a useful spectrum could be deblended for every bright star in the field of view. 



crease again, indicating that we have reached an "overcrowding 
regime" where our approach breaks down. 

We note that the distributions shown in Fig. [14] can also 
be very helpful to judge whether single-star spectra can be de- 
blended in a certain crowded stellar fields and to specify what 
spatial resolution would be required for its investigation. 

6. Potential sources for systematic errors 

6.1. Influence of the PSF 



In Sect. 4.2 we discussed to what accuracy the PSF can typically 
be determined in a crowded field datacube. We now investigate 
the effects of any possible mismatch between the true and the 
reconstructed PSF on the quality of the deblended spectra. To 
this aim, we modified the analysis of our two- star datacubes. The 
PSF assumed in the deblending process did now systematically 
deviate from the one that that was used to create the cubes. We 
describe the mismatch in terms of the FWHM of the PSF, as this 
is a quantity that is relat ively easily accessible in observations. 
The analysis of Sect. [42] showed that we can recover the FWHM 
to an accuracy of usually -1-2%, with some outliers at the 10% 
level. Based on these results, we simulated PSF mismatches of 
1,2, 5, and 10% in FWHM. 

Qualitatively, we expect the effect of an inaccurate PSF to 
be the following: The deblending process will leave residuals in 
the vicinity of every star. The amplitude of such residuals will 
scale with the brightness of the star. Close-by stars will then be 
deblended on top of the residuals and the extracted spectrum will 
be a combination of the true spectrum and the residuals. 

Quantitatively we wanted to know by how much these resid- 
uals bias a deblended spectrum and the derived astrophysical 
quantities. We measured again the radial velocities and Ca triplet 
equivalent widths in the deblended spectra, took the differences 
to the input values and checked the ensemble of results for sys- 
tematic deviations. 

For obvious reasons the impact of any PSF mismatch will 
be strongest for relatively faint sources in the vicinity of signif- 
icantly brighter ones. We therefore considered two cases: (i) A 
moderate brightness contrast between the source in question and 



its neighbour (1 < Am < 3), and (ii) a strong brightness con- 
trast of 3 < Am < 5. A contrast of 5 mag would be roughly the 
expected value for a star in the red clump of a Globular Cluster 
apparently close to a star at the tip of the red giant branch. 

In Fig. [15] we present the median absolute difference be- 
tween the measured and true radial velocities as a function of 
the degree of PSF mismatch, in two panels corresponding to the 
different contrast classes. To discriminate between random and 
possible systematic errors, the median absolute difference is also 
shown for the case case of a perfect PSF (red line in Fig. [15]). In 
this case, the offset should be completely caused by random er- 
rors, i.e. the limited S/N of the deblended spectra. Any increase 
in the offset can then be attributed to the imperfect PSF. With 
increasing influence of the PSF residuals of the brighter neigh- 
bouring star, the median of the distribution will approach the 
assumed velocity dispersion. 

As expected, the influence of the PSF increases with the 
brightness contrast between source and neighbour. This can be 



verified by comparing with the two panels of Fig. 15 For mod- 
erate brightness contrast, the introduced systematics are small if 
the FWHM of the PSF is determined to an accuracy of <5%, 
whereas in the case of a strong brightness contrast, offsets of 2% 
already introduce measurable systematics. At this contrast level, 
the residuals caused by the PSF mismatch if the errors in the 
FWHM are >5% are so strong that the signal of the fainter star 
basically disappears and no useful spectrum can be deblended 
any more. On the other hand, Fig. 15 also shows that PSF mis- 
match only becomes an issue for source separations comparable 
to or smaller than the FWHM. We find similar results regarding 
the accuracy of the recovered values of EWc a T- 

Recall that for crowded fields in Galactic globular clusters 
we typically can recover the PSF width to an accuracy of <2%. 
The simulations presented in this subsection demonstrate clearly 
that this will be sufficient to deblend an unbiased spectrum of 
a star in the close vicinity, to even a small fraction of the PSF 
width, of a neighbour that is ~10x brighter. Only in the extreme 
case where the brightness contrast between the two stars is sig- 
nificantly larger than a factor of 10, significant biases are to be 
expected. Yet such cases will be known from the input catalogue 
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Fig. 15. Difference between recovered and true radial velocity of a 
simulated apparent binary star, as a function of source separation and 
of the degree of PSF mismatch, for two ranges of flux ratios between 
the two stars. The curves show the median of the absolute velocity 
difference of each stellar pair. Also shown is the median offset for the 
case of a perfect PSF (red line) to show the behaviour expected in the 
absence of any systematic error. For clarity, deviations larger than the 
velocity dispersion have been set to a value of 10km s _1 , the assumed 
velocity dispersion in the cluster. 



and can thus be easily flagged and excluded from the further 
analysis. 

6.2. Further sources for systematic errors 

An imperfect PSF is not the only potential source of systematic 
errors. Another such source are the positions of the individual 
stars in the field. For reasons like measurement errors or proper 
motions, the source coordinates that are given in the input cata- 
logue can be offset from their true values. The strength of such 
effects is essentially a property of the specific input catalog that 
is used and cannot be predicted like the achieved accuracy of the 
PSF. For this reason, we did not try to quantify the influence of 
these effects on the deblended spectra. 

A further error source related to the input catalogue is the is- 
sue of missing stars or spurious detections. Those effects will oc- 
cur mainly in the regime of the fainter stars that are inaccessible 
to the integral-field observations (that have a lower spatial reso- 
lution). Yet under certain circumstances, they might also play a 
role amongst the brighter stars that are accessible. One such ex- 
ample are catalogues that are compiled from observations with 
the Advanced Camera for Surveys (ACS) onboard HST. In the 
case of globular clusters, those observations are usually targeted 
at the numerous faint main-sequence stars and the brightest gi- 




Fig. 16. Simulated MUSE observation of the globular cluster 47Tuc. 
Left: Cut-out from an HST/ACS observation of the central arcmin of 
the cluster in the F606W-passband. A red cross denotes the cluster 
centre. Right: Reconstructed broadband image from the mock MUSE 
data of the same region, obtained by integrating the datacube with the 
F606W filtercurve. The seeing in the simulation was set to 0.8arcsec. 
Red squares indicate the two 20 x 20 arcsec fields that are discussed in 
the text. 



ants appear heavily saturated in the exposures and cause strong 
bleeding features on the CCD that might cover relatively bright 
stars (see Fig [16] for an example). Again this effect will largely 
depend on the quality of the used input catalogue. 



7. MUSE 



MUSE ( [Bacon et al.|2010| ) is an integral field spectrograph cur- 
rently being built by a consortium of 6 European institutes and 
ESO. It is scheduled to see its first light at the Very Large 
Telescope (VLT) in 2013. The instrument provides a FoV of 
1 arcmin 2 with spaxels of 0.2 x 0.2 arcsec 2 , and a wavelength 
coverage of 4650-9300A. The combination of a large FoV with 
a spatial sampling sufficient to properly sample the PSF even 
under good seeing conditions makes MUSE a unique instrument 
for a variety of science applications. Although the main moti- 
vation for developing this new instrument is the observation of 
faint galaxies at medium to high redshift, some very promising 
applications exist for the investigation of crowded stellar fields. 
To demonstrate this, we outline in the following the analysis of 
a simulated MUSE datacube of the Globular Cluster 47Tuc. 

The simulations are again based on the HST-photometry ob- 
tained in the HST/ACS Survey of Galactic Globular Clusters. 
Based on broadband colours and an isochrone fit to the colour 
magnitude diagram of 47Tuc, each star was assigned a spectrum 
based on a new library of model atmospheres and synthetic spec- 
tra calculated by Husser et al. (submitted) using the stellar atmo- 
sphere code PHOENIX ( [Hauschildt & Baron|1999] ). 

To simulate the effect of missing stars in the vicinity of 
brighter stars, we applied the following correction to the input 
catalogue: we counted the surface density of stars at a given 
magnitude in the vicinity of brighter stars and compared it to the 
overall density of those stars across the field covered by the cat- 
alogue. Stars were then randomly added in the vicinity of the 
brighter ones until the two densities matched.The catalogue that 
was later used in the analysis did not include those stars. 

In the final step of the simulation, a datacube was created us- 
ing dedicated software developed within the MUSE consortium 
(R. Bacon, private communication). It creates a datacube con- 
taining the provided sources and a sky spectrum. Each spectrum 
is convolved with the line spread function of MUSE. The seeing 
in the simulation was set to 0.8 arcsec. This value is internally 
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Fig. 17. Visualization of our source deblending algorithm applied to simulated MUSE data. For each of the two fields highlighted in Fig. 16 
we show a cut-out from an HST/ACS image (left), a white light image of the simulated data (centre), and a residual image after fitting the sources 
(right). Each location where a useful spectrum was deblended has been marked by a blue circle. 



translated to a wavelength dependent PSF. The final simulated 
datacube is a combination of three snapshot exposures, each with 
an exposure time of 30s. Fig.[l6]shows a whitelight image of this 
datacube together with an HST image. 

We concentrate our discussion on two 20 x 20 arcsec sub- 
cubes, highlighted in Fig. [16] (right) by red squares. In both re- 
gions, we deblended the stellar spectra using our crowded field 
spectroscopy code. 

The algorithm itself was slightly modified in order to handle 
the significantly larger amount of data of a MUSE cube com- 
pared to PMAS or ARGUS. The size of the FoV of MUSE is suf- 
ficiently large so that some relatively isolated bright stars should 
exist within the FoV that can be used as PSF calibrators. We 
therefore used the existing photometry to select suitable PSF 
stars on the condition that within a given radius there are no 
neighbours brighter than a given flux ratio. The search radius 
is chosen such that the PSF contribution outside this radius is 
essentially zero. After each fit to the object fluxes during the 
iteration, we then subtract all stars except those previously iden- 
tified and determine the coordinate transformation and the PSF 
using only those stars. This significantly speeds up the analy- 
sis and we achieve computational times similar to our PMAS or 
ARGUS datacubes. Nevertheless, the actual extraction of the 
spectra still is performed on all stars simultaneously. 

In Fig. [T7] we show a close-up of the two regions, again 
using an HST/ACS image and a whitelight image of the mock 
MUSE data. We also show a whitelight image for each region 
with the MUSE data where the deblended sources were already 



subtracted. Closer inspection of these residuals reveals that some 
stars have been missed by our source selection; these are the stars 
added to the incomplete source catalogues as discussed above. 
Such sources can be easily identified in the residuals and then 
added manually to the catalogue. 

To visualize the efficiency of our deblending approach we 
marked the position of every source for which a useful spectrum 
was deblended. Note that sources for which the extracted spectra 
have a S/N too low for a reliable radial velocity determination 

The total number of useful spectra 



17 



are not marked in Fig. 
that were deblended is 580 in subfield #1 and 610 in subfield #2. 
Interpolating these numbers to a full datacube, we estimate that 
from a single MUSE observation obtained under average seeing 
conditions we can obtain ~ 5000 useful spectra; under very good 
seeing conditions this number may even be 3-4x higher. 

The extended wavelength range of MUSE allows us to di- 
rectly determine broadband colours from the spectra by apply- 
ing the filtercurves of the HST/ACS F606W and FUAW filters 
(hereafter called V and /, respectively) and compare them to the 
'true' (i.e., input) colours. In Fig.[l8j the deviation in V - I is 
plotted as a function of V band magnitude. The comparison of 
the two fields gives a good impression of the effect of crowd- 
ing: In the less crowded part of the simulated data (field #1), we 
obtain useful spectra for fainter stars than in the direct vicinity 
of the cluster centre (field #2). Yet in both fields we are able 
to probe below the main sequence turn-off. Taking into account 
that the simulation assumed only average seeing conditions and 
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Fig. 18. Deviation in V-I colour recovered from the deblended spectra 
as a function of a stellar magnitudes for the two subfields highlighted 
in Fig. [16] The median of the distribution is shown as a thick solid line, 
dashed lines give the 75% percentiles of the distribution. A vertical 
dashed line is used to indicate the main-sequence turn-off in both panels. 



Fig. 19. Deviation of the radial velocities determined from the de- 
blended spectra, again plotted as a function of a stars brightness for the 
two fields highlighted in Fig. [16] Radial velocities were determined by 
cross-correlating each deblended spectrum with its input spectrum. The 
line types are as in Fig. [18] 



that MUSE is also designed to work with adaptive optics, this 
demonstrates the unique capabilities provided by the instrument. 

Finally, we take a look at the accuracy achievable in the 
measured radial velocities. The spectral resolution provided by 
MUSE is smaller than what we used in our previous simulations 
and that of our existing PMAS and ARGUS data. However, the 
larger wavelength range at least partly compensates for this. As 
Fig. 19] shows, radial velocities of bright giants can be deter- 
mined with an accuracy of a few km/s. Around the main se- 
quence turn-off of the cluster, the typical error of the obtained 
radial velocities is still comparable to the velocity dispersion. 

8. Conclusions 

The application of PSF fitting techniques to integral field spec- 
troscopy is a powerful approach to observe crowded stellar 
fields. We developed an algorithm to deblend the spectra of 
many stars within a single datacube simultaneously, and vali- 
dated the method by applying it to realistic simulated observa- 
tions of the central regions of globular clusters. The combination 
of linear least- squares fitting with the usage of sparse matrices 
makes the code computationally efficient and affordable even for 
a modest workstation. 

One central assumption for our algorithm is that an input cat- 
alogue of the sources in the field already exists. Typically this 
catalogue would be obtained by other means, such as by running 
a classical crowded field photometry code on high-resolution 
HST images. It is of course also possible to perform the source 
detection in the IFS datacubes themselves, e.g. from a collapsed 
white-light image. But our simulations show that by using prior 
knowledge of the locations of sources in the field, the number of 
correctly deblended sources is increased by up to a factor of ~4 
compared to the case where the source detection is performed on 
the IFU data. 



We have extensively tested the performance of our code as 
a function of the degree of crowding, expressed as the number 
of sources per field or per spatial resolution element. A con- 
ventional rule of thumb for crowded field photometry states that 
deblending performs well up to a stellar density of ~0. 1 per res- 
olution element. We have shown that the spectroscopic deblend- 
ing works at even considerably higher source densities than that. 
This gain is partly due to the application of prior knowledge as 
discussed above, and partly due to the continuity enforcement 
over many simultaneously evaluated image layers in a datacube. 

While unbiased spectra can be extracted even for heavily 
blended sources, such spectra will suffer from a significantly re- 
duced S/N level, with the degradation being driven by the prox- 
imity to and the brightnesses of nearby stars. We have shown 
that this reduction of S/N due to blending can be accurately mod- 
elled and predicted for a given dataset from the input catalogue. 
Consequently, an optimal source list for the final extraction can 
be constructed according to the expected S/N of the final spec- 
tra. This is a very useful feature for statistical investigations in 
crowded stellar fields, as it allows one to maximise the number of 
'meaningful' spectra that can be obtained from a given dataset. 

Under conditions of strong crowding, the number of stars per 
field for which spectral parameters can be reliably determined 
is approximately independent of the actual source density and 
corresponds to roughly 0.2 stars per spatial resolution element. 
This 'plateau' exists because of the mutually opposing effects of 
higher source densities on the one hand, and more severe S/N 
degradation due to crowding on the other hand. Only when the 
observed density of stars of comparable brightness passes a def- 
inite 'overcrowding limit' of ~1 star per resolution element, the 
extraction of useful individual spectra breaks down entirely. 

The degree of crowding in a given field depends of course 
also on the depth and angular resolution of the data. We con- 
structed our simulated datacubes in view of our own existing ob- 
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servations of Galactic globular clusters, using present-day IFUs 
under seeing-limited conditions. Even in the very central regions 
of these clusters we found that the source density of stars bright 
enough to produce meaningful spectra was still well below the 
overcrowding limit. However, this could change rapidly if the 
data were going deeper down the colour-magnitude diagram, 
especially once the main sequence is being probed. We have 
demonstrated that with the upcoming MUSE instrument this do- 
main will actually be reached. Improving the angular resolution, 
for example through ground-layer adaptive optics as envisaged 
for MUSE, will then become crucial. 

The spectroscopy of crowded stellar fields may be of inter- 
est for other classes of astronomical objects, such as compact 
open clusters, dwarf galaxies, or dense regions in the bulge of 
the Milky Way. The methodical work presented in this paper 
will enhance the capabilities of 'crowded field 3D spectroscopy' 
beyond our own application topic of globular clusters. For the 
benefit of the community, we plan to make our code available to 
the public in the future. 
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